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1  Introduction 

The  overall  goal  of  our  AFOSR  sponsored  research  program  is  to  construct  fast  and 
efficient  numerical  algorithms  for  solving  stochastic  partial  differential  equations  and 
apply  them  to  solve  optimal  control  problems  under  uncertainty  which  is  described 
by  stochastic  partial  differential  equations. 

Uncertainty  quantification  has  been  an  active  research  area  in  the  past  10  years 
with  many  significant  applications  ranging  from  signal  processing  to  aircraft  wing 
designs.  It  is  well  understood  that  effective  numerical  methods  for  stochastic  partial 
differential  equations  (SPDES)  are  essential  for  uncertainty  quantification.  In  the 
last  decade  much  progress  has  been  made  in  the  construction  of  numerical  algorithms 
to  efficiently  solve  SPDES  with  random  coefficients  and  white  noise  perturbations. 
However,  high  dimensionality  and  low  regularity  are  still  the  bottleneck  in  solving  real 
world  applicable  SPDES  with  efficient  numerical  methods.  This  project  is  intended  to 
address  the  mathematical  aspects  of  numerical  approximations  of  SPDES,  including 
error  analysis  and  complexity  analysis  and  development  of  new  efficient  numerical 
algorithms.  Our  contributions  toward  this  objective  include 

(i).  Optimal  Kronecker  sequence  for  the  quasi  Monte  Carlo  method  and  the  Monte 
Carlo  method  for  stochastic  Stokes  equations.  The  generation  of  appropriate  high- 


quality  quasirandom  sequences  (low-discrepancy  sequences)  is  crucial  to  the  success 
of  quasi-Monte  Carlo  methods.  We  present  a  new  algorithm  for  hireling  an  optimal 
Kronecker  sequence  within  choices  of  irrationals,  we  illustrated  with  numerical  ex¬ 
periments  why  onr  algorithm  is  efficient  for  breaking  these  correlations. 

(ii) .  Efficient  algorithms  for  stochastic  partial  differential  equations.  We  are  inter¬ 
ested  in  the  the  orthogonal  polynomial  expansion  of  d-variable  functions.  Fast  and 
efficient  algorithms  of  evaluating  such  expansions  are  essential  for  efficiently  solving 
stochastic  partial  differential  equations  with  spectral  methods  where  d  may  be  fairly 
large. 

(iii) .  High  order  numerical  method  for  nonlinear  filter  problems.  We  first  convert 
the  nonlinear  filter  problem  into  a  problem  of  solving  a  stochastic  partial  differential 
equations  (SPDES).  We  then  construct  a  high  order  method  to  solve  the  SPDES  and 
apply  our  algorithm  to  solve  a  ship  tracking  problem. 

2  Brief  Overview  of  Accomplishments 

During  the  prior  three-year  grant  period  we  have  published  12  papers  that  have 
been  supported  by  AFOSR  funding.  Selections  of  onr  current  accomplishments  are 
summarized  here  to  provide  an  introduction  to  the  concepts  we  have  developed  in  the 
past  and  upon  which  future  planned  studies  are  based. 

2.1  The  optimal  Kronecker  Sequence  for  quasi  Monte  Carlo 
simulations 

Kronecker  sequences  are  easy  and  fast  to  implement.  However,  problems  with  Kro¬ 
necker  sequences  arise  from  correlations  between  the  choice  of  parameters  which  are 
used  for  different  dimensions.  These  correlations  cause  Kronecker  sequences  to  have 
poor  two-dimensional  projections  for  some  pairs  of  dimensions. 


Figure  1:  Left:  1000  points  of  dimension  1  and  2  projection  for  original  Kronecker 
sequence;  right:  1000  points  of  dimension  26  and  27  projection  for  original  Kronecker 
sequence. 

In  Figure  1,  we  can  see  high  correlations  and  similarities  in  the  poor  two-dimensional 
projections  for  original  Kronecker  sequences. 

There  are  at  least  two  possible  ways  to  break  the  correlations  we  have  seen  in 
the  Kronecker  sequence.  One  is  by  increasing  the  difference  between  the  bases  for 
any  pair  of  dimensions;  the  other  is  to  select  better  basis  for  the  Kronecker  sequence. 
It  became  clear  that  the  smaller  the  partial  quatients  in  the  coninued  fractions  of 
the  irrarional  number,  the  more  uniformly  distributed  the  s-dimcnional  Kronecker 
sequence  is  [4],  We  blend  those  two  methods  and  propose  a  novel  approach  by  using 
generarlized  golden  ratio,  which  is  widely  applied  in  many  applications  in  computer 
science  [6].  The  2-D  projections  derived  from  our  new  algorithm  are  shown  in  Figure 
2.  The  patterns  appeared  in  Figure  1  disappear  in  Figure  2. 

We  have  demonstrated  that  generalization  of  golden  ratio  is  the  simplest  and  most 
effective  method  to  break  this  correlation  in  original  Kronecker  sequences.  It  is  also 
very  easy  to  implement  the  original  Kronecker  sequences  and  generate  high-quality 
of  the  optimal  Kronecker  sequences. 

In  Bratley  et  al’s  paper  [1],  the  Faure  and  Sobol  sequences  are  used  to  evaluate  high 
dimensional  integrals,  and  the  errors  in  the  numerical  results  for  over  30  dimensions 
become  quite  large.  When  we  compared  their  results  with  our  optimal  Kronecker 
sequences,  we  found  that  our  method  has  much  smaller  errors  comparing  to  the 


Figure  2:  Left:  1000  points  of  dimension  1  and  2  projection  for  optimal  Kronecker 
sequence;  right:  1000  points  of  dimension  26  and  27  projection  for  optimal  Kronecker 
sequence. 

benchmark  Faure  and  Sobol  methods. 

2.2  High  order  numerical  method  for  nonlinear  filter  prob¬ 


lems 


Assume  that  X  =  Xt  is  a  signal  process  and  Y  =  Yt  is  the  observation.  Because  of 
noises,  X  and  Y  are  solutions  of  the  following  stochastic  differential  equations. 


dXt  =  b{Xt)dt  +  a{Xt)dUt 
dYt  =  h(Xt)dt  +  dVt 


ffere  Ut,  Vj  are  independent  Brownian  motions  which  reflect  the  noises  of  the  signal 
and  the  observation.  The  Liter  problem  is  to  seek  the  best  estimation  of  (j>{X)  given 
the  observation  Y.  Mathematically  this  amounts  to  seek  a  conditional  expectation 


E  (<f>(Xt)  |  (Ys,  0  <  s  <  t)),  i.e., 


E  (<f>(Xt)  |  (ys,0  <s<t))  =  Argmin  (E  \</>(Xt)  -  Zt\2) 

{Zt  is  Y(Q  t)  measurable} 


According  to  Zaikai  Liter  theory,  the  best  approximation  is  given  by 


E  ((/)(Xt)  | 


Figure  3:  Sample  path  of  the  ship  location 


Figure  4:  Un-normalized  probability  distribution  of  the  ship  location 


Also  according  to  the  Zakai’s  theory,  p  satisfies  the  following  stochastic  partial  differ¬ 
ential  equation. 


dpt  =  {  aa*  PJ  +  V]  b~-)dt  +  pth*dYt. 

OXiOXj  r)'r 


i,j= 1 


i= 1 


So  the  nonlinear  filter  problem  is  solved  as  long  as  the  above  stochastic  partial  dif¬ 
ferential  equation  is  be  solved.  We  have  derived  a  higher  order  method  to  solve 
the  nonlinear  stochastic  differential  equations.  We  achieve  the  high  order  method 
through  solving  a  backward  forward  doubly  stochastic  differential  equations.  Detail 
of  the  method  can  be  found  in  [3]. 

We  then  applied  the  method  to  solve  a  problem  of  tracking  a  ship  in  open  water 


through  the  angle  observation  of  the  ship.  Figure  3  shows  a  sample  path  of  the  ship 
and  Figure  4  shows  the  probably  distributions  of  the  ship  location  at  four  different 
times. 


2.3  Error  analysis  of  finite  element  approximations  for  stochas¬ 
tic  Stokes  equations 

In  this  subproject,  we  consider  the  steady-state  Stokes  equation  with  the  forcing  term 
perturbed  by  white  noise: 

— uAu(x,  £)  +  Vp(x,  £)  =  f  +  cr(x)W(x,£)  x  G  DcRd(d  =  2,3),  £  e  fl, 

<  divw  =  0  in  D, 

u  —  0,  on  d D. 

(1) 

Here  (fi,  E,  P)  is  a  probability  space  and  W  =  (W1, . . . ,  Wd)  is  the  white  noise  such 
that 

E{Wi(x)Wi{x?))  =  5{x  -x')r  x,  x'  e  fi,  j  =  1, . . . ,  d. 

Our  strategy  of  solving  this  problem  consists  of  two  steps.  In  the  first  step,  we 
construct  an  approximate  white  noise  Ily,  based  the  finite  element  partition  of  the 
physical  domain  D.  In  the  second  step,  we  construct  a  finite  element  approximation 
(uhiPh )  f°r  the  solution  (u,p)  of  (1).  Under  certain  assumptions  on  the  regularity  of 
the  domain  and  using  the  Taylor-Hood  finite  element,  we  proved  the  following  error 
estimates  [2], 

f  E(\\u-uh\\2+  Wp-PhW2^)  <  C\\nh\h2  d  =  2, 

|  E  (||w  -  uh ||2  +  ||p  -  PhW-i)  <Ch  d  =  3. 

The  significance  of  the  error  estimate  is  twofold: 

•  It  indicates  that  the  due  to  the  lack  of  regularity,  the  convergence  order  of  the 
finite  element  solution  for  stochastic  Stokes  equations  is  much  lower  than  that 
for  deterministic  Stokes  equations. 

•  It  provides  a  guidance  for  the  number  of  samples  used  in  Monte  Carlo  simula¬ 
tions  to  evaluate  the  statistics  of  the  solutions. 


We  also  constructed  practical  numerical  algorithms  to  find  the  finite  element  approx¬ 
imations  using  the  Monte  Carlo  method  and  verified  our  theoretical  results  through 


numerical  experiments  ([2]). 


2.4  Fast  and  efficient  numerical  algorithms  for  stochastic  par¬ 
tial  differential  equations 

We  first  developed  a  fast  algorithm  of  high  dimension  orthogonal  polynomial  ex¬ 
pansions  on  sparse  grids.  Then  we  apply  this  algorithm  to  solve  stochastic  partial 
differential  equations  (SPDEs).  This  is  achieved  through  the  following  four  steps. 

•  Step  1.  Use  only  sparse  indices  in  the  orthogonal  expansions.  For  a  d  dimen¬ 
sional  problem,  this  reduces  the  number  of  terms  in  the  orthogonal  expansion 
from  0(nd )  to  0(nlogd_1n),  a  tremendous  saving  of  computing  cost. 

•  Step  2.  Use  fast  Fourier  transform  as  well  as  sparse  grid  to  compute  the  coeffi¬ 
cients  of  the  orthogonal  expansion.  This  reduces  the  computing  cost  from  0(nd) 
to  0(n  log(d+1)n\  yet  another  significant  reduction  of  computing  complexity. 

•  Step  3.  Apply  the  fast  orthogonal  expansion  algorithm  to  solve  stochastic 
partial  differential  equations  with  random  coefficients  and  forcing  term  using  the 
spectral  method.  The  remarkable  feature  of  our  algorithm  is  that  the  coefficient 
matrix  is  sparse  and  each  element  can  be  calculated  analytically. 

Our  numerical  experiments  show  that  the  overall  computing  complexity  of  solving 
the  stochastic  partial  differential  equation  is  in  the  same  magnitude  as  computing 
the  orthogonal  expansion  of  the  forcing  term,  which  is  proved  to  be  0(n\ogd+l  n). 
In  comparison,  the  complexity  of  conventional  spectral  method  is  at  least  0(nd), 
which  is  prohibitively  expensive  when  the  dimensional  d  is  large  as  in  the  case  of 
solving  stochastic  partial  differential  equations.  The  numerical  verification  of  the 
computational  complexity  is  shown  in  Figure  2.4. 


Figure  5:  The  CPU  times  in  numerical  experiment  v.s.  theoretical  estimates  Solid 
lines:  0(nlogd+1  n),d  =  4,6,8  Dashed  lines:  CPU  times  in  numerical  experiment 
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